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^ , Abstract 

\^ ' This paper collects the efltorts done in our previous works [8].|llj.[T0] to build a robust 

, multiscale kinetic-fluid solver. Our scope is to efficiently solve fluid dynamic problems which 

present non equilibrium localized regions that can move, merge, appear or disappear in time. 
I The main ingredients of the present work are the foUowings ones: a fluid model is solved in the 

QQ . whole domain together with a localized kinetic upscaling term that corrects the fluid model 

' wherever it is necessary; this multiscale description of the flow is obtained by using a micro- 

^\ , macro decomposition of the distribution function [10] : the dynamic transition between fluid 

' and kinetic descriptions is obtained by using a time and space dependent transition function; 

to efficiently define the breakdown conditions of fiuid models we propose a new criterion based 
on the distribution function itself. Several numerical examples are presented to validate the 
' method and measure its computational efficiency. 

' 

■ Keywords: kinetic-fluid coupling, multiscale problems, Boltzmann-BGK equation. 



1 Introduction 

Many engineering problems involve fluids in transitional regimes such as hypersonic flows or micro- 
electro-mechanical devices. In these cases, usual fluid models (like Euler or Navier-Stokes equa- 
tions) break down in localized regions of the computational domain (typically in shock and bound- 
ary layers). For such problems, using classical fluid models is generally not sufficient for an accurate 
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description of the flow in non-equilibrium regions. However, it is not necessary to solve the Boltz- 
mann equation — which is computationally more expensive than continuum solvers by several orders 
of magnitude — especially in situations where the flow is close to thermodynamical equilibrium. 

For the above reasons, it is important to develop hybrid techniques which can reduce the 
use of kinetic solvers to the regions where they are strictly necessary (kinetic regions) , leaving the 
simulation in the rest of the domain to a continuum or fluid solver (fluid regions). The construction 
of these methods involve two main problems. The flrst one is how to accurately identify the different 
regions. We refer, for instance, to the works of Wijesinghe and Hadjiconstantinou [15 , Levermore, 
Morokoff, and Nadiga [20,, and Wang and Boyd 32J, in which various breakdown criteria are 
proposed. The second main problem is how to efficiently and correctly match the two models 
at the interfaces. Most of the recent methods are based on domain decomposition techniques, 
such as in the works of Bourgat, LeTallec, Perthame, and Qiu [3j, Bourgat, LeTallec and Tidriri 

LeTallec and Mallinger [53], Aktas and Aluru [T], Roveda, Goldstein and Varghese ^7\, Sun, 
Boyd and Candler ^25], Wadsworth and Erwin f33', and Wijesinghe et al. [Mj. The same domain 
decomposition approach has been also used in many others fields, such as, for instance, in molecular 
dynamics [M], in epitaxial growth [28] or for problems involving diffusive scalings [18] instead of 
hydrodynamic ones. We also mention the use of decompositions in velocity instead of physical 
space done by Crouseilles, Degond and Lemon \Z, and by Dimarco and Pareschi |12j . 

It is important to stress that most of the mentioned methods use a static interface between 
kinetic and fluid regions that is chosen once for all at the beginning of the computation. However, 
for unsteady problems, this approach appears as somehow inadequate and inefficient, and for this 
reason, some automatic domain decomposition methods have also been proposed, see for example 
Kolobov et al. [^ , or Tiwari [301 131] and Dimarco and Pareschi fTP . We have also proposed a 
similar approach in 

In this paper, we propose a method that has similar features as the methods mentioned above: 
we solve the Boltzmann-BGK equation coupled with the compressible Euler equations through an 
adaptive domain decomposition technique. With this technique, it is possible to achieve consid- 
erable computational speedup, as compared to steady interface coupling strategies, without losing 
accuracy in the solution. This method is somehow an extension of our previous work [llj . but 
several important differences must be noted. First, we introduce a new breakdown criterion which 
is based on a careful inspection of the distribution function. This criterion can be defined by using 
the macroscopic variables only, at least in fluid regions, and thus does not introduce additional 
expensive computations. This allows us to define kinetic regions that are small as possible. Second, 
we use a decomposition of the distribution function that has better properties than the one used 
in [11]. In fact, while it has been proved by Degond, Jin, and Mieussens [8] that the decomposition 
used in [11] preserves uniform flows at the continuous level, we show in this paper that this is not 
true in general at the discrete level, except if a quite specific and very expensive kinetic scheme 
is used. For this reason in the present work we use the decomposition proposed by Degond, Liu, 
and Mieussens [lOj . since it perfectly preserves uniform flows, both for the continuous and discrete 
cases. As in [10], we decompose the distribution function into an equilibrium part, that can be 
described by macroscopic fluid variables, and a perturbative non-equilibrium part. We obtain a 
micro-macro fiuid model in which the macroscopic variables are determined by solving a fiuid equa- 
tion with a kinetic upscaling. This kinetic upscaling is determined by solving a kinetic equation, 
and is dynamically and automatically localized wherever it is necessary, by using our breakdown 
criterion and the transition function idea [Q] [H [TOl [11] . Third, we propose an efficient numerical 
scheme for discretizing our micro-macro fluid model: we use a time splitting approach that has 
several advantages. In particular it is shown to preserve the positivity of the distribution function. 

The outline of the article is the following. In section [2l we introduce the BGK equation and 
its properties. In section jS] we present the coupling strategy, while in section [4] the numerical 
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scheme is described and positivity properties are analyzed. In section [SI we derive our breakdown 
criterion, and the final algorithm is presented. Several numerical tests are presented in section [5] 
to illustrate the properties of our method and to demonstrate its efficiency. A short conclusion is 
given in section [71 In appendix A, the differences between the present coupling strategy and the 
decomposition used in [8l lllj are analyzed in some details. 



2 The Boltzmann-BGK model 

We consider the kinetic equation 

dtf + vy.f = Qif), (1) 

with the initial data 

f{x,v,t^ 0) = /™t, 

where / — f{x, v, t) is the density of particles that have velocity u e and position x G C M'^ 
at time t > 0. The collision operator Q locally acts in space and time and takes into account 
interactions between particles. It is assumed to satisfy local conservation properties 

(mQ(/)) = (2) 

for every /, where we denote weighted integrals of / over the velocity space by 

m = / Hv)f{v)dv, (3) 



where (f>{v) is any function of v, and m{v) = {l,v, \v\'^) are the so-called collisional invariants. It 
follows that the multiplication of ([T|) by m{v) and the integration in velocity space leads to the 
system of local conservation laws 

at(m/) + V,(wm/) =0. (4) 

We also assume that the functions satisfying Q{f) = 0, referred to as local equilibrium distributions 
and denoted by E[q\, are defined implicitly through their moments g by 

Q = {mE[Q]) (5) 

In the present paper we will work with the BGK model of the Boltzmann collision operator that 
reads 

Qif)^iyiE[e]-f). (6) 

With this operator, collisions are modelled by a relaxation towards the local Maxwellian equilib- 
rium: 

g / — |w — 



^[^](^)-(^-p^^^J, (7) 

where g and u are the density and mean velocity while 9 = RT with T the temperature of the gas 
and R the gas constant. The macroscopic values g, u and T are related to / by: 



Q - 



I fdv, gu^ ( vfdv, = ^1 \v- u\^fdv, (8) 



while the internal energy e is defined as 



e=^ / \vffdv^l\uf + h. (9) 
Zg Jr3 Z Z 
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The parameter i/ > is the relaxation frequency. In this paper, we use the classical choice v = /i/p 
where \i = ^ref{(^/Oref)^ is the viscosity and p is the pressure. We refer to section [B] for numerical 

values of /^re/, dref and LjJ. 

Boundary conditions have to be specified for equation ([T]). Different type of conditions are 
used in applications: inflow, outflow, specular reflection or total accomodation. We will specify 
the conditions we use for every numerical test in section [6l 

When the mean free path between particles is very small compared to the size of the compu- 
tational domain, the space and time variables can be rescaled to 

x' = ex, t' = et (10) 

where e is the ratio between the microscopic and the macroscopic scale (the so-called Knudsen 
number). Using these new variables in ([T]), we get 

dt'r + vy.'r = '^iE^Q]-n. (n) 

If the Knudsen number e tends to zero, this equation shows that the distribution function converges 
towards the local Maxwellian equilibrium E^[q\. Using this relation into the conservation laws ^ 
gives the Euler equations for g: 

dt'Q + V.'F{Q)^Q, (12) 

where F{q) — {vmE^[Q\) — {g, gu, ge). 

In the sequel, to give a simple description of our approach, all schemes and all algorithms 
are shown for the one dimensional case in velocity and physical space. The extension to the 
multidimensional case does not introduce any additional difficulty in the mathematical setting. 
We will also omit the primes wherever they are unnecessary. 

3 The coupling method 

In this section, we follow the work of [10] and extend the micro-macro fluid model to allow for 
dynamic localization of the kinetic upscaling. 

3.1 Decomposition of the kinetic equation 

Our method is based on the micro-macro decomposition of the distribution function: it is decom- 
posed in its local Maxwellian equilibrium and the deviation part as 

f = E[Q]+g. (13) 

Because the equilibrium distribution has the same first three moments as / we have 

{mg) = 0. (14) 

Then it can be easily proved that the following coupled system 

dtQ + d^F{Q)+d^{vmg) ^0 (15) 
dtg + vdxg ^ -vg ~ {dt + vdx)E[Q] (16) 

is satisfied by ^ = (™/) and g = f — E[q], where F{g) = {vmE[Q]) is the flux associated to the 
equilibrium state. The corresponding initial data are 

Qt=() = Qtmt = {■mf^n^t), 9t=Q = finit ~ ^^iffimt]- (17) 
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The converse statement is also true: if q and g satisfy system and with initial data (|17p. 
then / = E[q\ -\- g satisfies the kinetic equation ([1]) (see [S] for details). In the following sec- 
tion, starting from this decomposition, we introduce the set of equations that define the domain 
decomposition technique we are proposing. 

3.2 Transition function 

Let fii, 0.2, and Q.^ be three disjointed sets such that fii U U = R'^. The first set fJi is 
supposed to be a domain in which the flow is far from the equilibrium (the "kinetic zone"), while 
the flow is supposed to be close to the equilibrium in VI2 (the "fluid zone") and also in (the 
"buffer zone"). We define a function h{x,t) such that 

{1, for X g rii, 

0, for X € rj2, (18) 

0<h{x,t)<l, for xefls. 

Note that the time dependence of h means that we account for possibly dynamically changing fluid 
and kinetic zones. The topology and geometry of these zones is directly encoded in h and may 
change dynamically as well. 

Next, we split the perturbation term in two distribution functions gK — hg and gp — {1 — h)g. 
The time derivatives of these functions then are 

dtgK = dt{hg) ^ gdth + hdtg, 

dtgp = dt{{l - h)g) = -gdth + (1 - h)dtg, 

and it is therefore easy to derive the following coupled system of equations 

dtQ + dxF{Q) + dx{vmgK) + dxivrngp) = (19) 

dtgK + hvdxgK + hvO^gp = -i^9k - h{dt + vdx)E[Q\ + ^^th, (20) 

QF 

dtgp + (1 - h)vd^gK + (1 - h)vd^gF = -vgp - (1 - h){dt + vd^)E[Q] - ^—j^dth, (21) 
with initial data 

gK,t=o = ht=ogt=o , 5^4=0 = (1 - ht=o)gt=o , et=o = Qinit (22) 

and with ht=o = hinit and gt=o = finit — -^'[^mit]- Again, system (jl9H2ip with initial data (j22p is 
equivalent to system (|15f[T6t with initial data (fT7|) (see [8] for details). 

Now assume that the flow is very close to equilibrium in f22 U fia. This means that g is very 
small in these domains and can be set to zero. Since g = gp in ^1,2, we set gp = in this domain. 
In 03, we also set gp — 0, which means that we approximate g by gK- Consequently, gp can be 
eliminated from (jT^HUl) to get: 

dtQ + d^F{Q) + d^{vmgK) = (23) 

dtgK + hvd^gK = --gK- h{dt + vd,)E[Q\ + ^dth, (24) 

with initial data 

gK,t=0 = ht=ogt=0 ^ hinitifmit ^ E[Q„^if]) , Qt=0 = Qinif (25) 

Note that since by definition gK is zero in the fluid zone ^2, the kinetic equation equation ([24]) is 
solved in the kinetic and buffer zones Sli and U,^ only. Indeed, in the fluid zone, we only solve (|23p 
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with gx — 0, which is nothing but the Euler equations. In the kinetic zone, we have — g 
and hence system (|23f[Ml) is nothing but system (fTBHTO)) . which is equivalent to the original BGK 
equation. System (|23H24p is our micro-macro fluid model with dynamically localized kinetic effects 
which will be used to solve multiscale kinetic problems. With this system, the distribution function 
/ is approximated by E[q\ + g^- 

In the next section we describe and analyze the numerical scheme we use to discretize this 
system, and we compare this new model to the model used in our previous work jllj . 

Remark 1. We mention here a slightly different derivation that leads to a different micro-macro 
model. In i20\} . the term ^ can be equivalently replaced by gx + gp, since gx — hg by definition 
and also g = gx + gp- In this case, the approximation gp = in 5^2 cind fia leads to the model: 

dtQ + d^F{Q) + d^vrngx) = (26) 

dtgK + hvd^gx = --gK - h{dt + vdx)E[Q] + gxdth (27) 

e 

Note that, surprisingly, this model is different from i2S[]24^ : indeed, the factor of dth is gx in h26[ - 
\27^ , while it is ^ in i23\\24\ l- However, we only use system i23\\24^ in the sequel, since it can be 
proved to have very good properties (like positivity preservation) . 

4 Numerical approximation of the coupled model and its 
properties 

First, we briefly describe a velocity discretization of the kinetic BGK equation. Then, we propose 
a second order in space numerical scheme for the perturbation term gx and for the macroscopic 
fluid equations. Then, we introduce a time splitting method between the transition function term 
and the rest of the system to compute the evolution of the perturbation function gx- Finally, 
positivity property for the distribution function / is analyzed in details. 

4.1 Discrete velocity model for kinetic equations 

Here, we replace the continuous velocity space by a bounded Cartesian grid V of iV nodes Vj = 
jAv + a, where j is a bounded index, Av is the grid step, and a is a constant. The collisional 
invariants m{v) are replaced by mj — {l,Vj, ^|wjP). The distribution function / is approximated 
on the grid by {fj{t,x))j, where fj{t,x) w f{x,Vj,t), while the fluid quantities are obtained from 
fj through discrete summations on V: 

3 

The BGK model is then replaced by the following system of N hyperbolic equations with a stiff 
relaxation term: 

= ^(£,[e]-/,), (29) 

where £j[Q] is the approximation of the continuous Maxwellian E[g\. Note that this approximation 
is not the evaluation of E[q\ on the grid points: in fact, to ensure conservation of macroscopic 
quantities and entropy decay at the discrete level, the approximated Maxwellian £[Qj] is defined 
through an entropy minimization problem that can be solved by computing the solution of a small 
non-linear system (we refer the reader to [5T1 [52] for details about this approximation) . 
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Finally, a micro-macro system with localized upscaling can be derived from the discrete- velocity 
BGK equation (|29p . exactly as in the continous case, and we find: 



dtQ + d^F{Q) + d^{vmgK) = (30) 
dtgK,] + hvjdxOKj = ~~9K,j - h{dt + Vjdx)£j[Q\ 4- -^dth, (31) 

with initial data 



9K.J,t=0 — ht=ogt=O.J — hinitifinitj - £][Qinit\) > £"(=0 " Qtnit' 

where (.) now stands for .Aw. 

4.2 Numerical schemes 
4.2.1 Non-splitting scheme 

For the space discretization, we consider a grid of step Ax and nodes Xi, while for the time 
discretization, we consider the step At and times t„ = nAt. The unknowns g and g^j are 
approximated by ~ g(tn,Xi) and gxij ~ gj{tn,Xi). Now, the space and time discretization of 
the discrete velocity micro-macro system (I30l - |3T|) is: 



(32) 



At ' \ Ax J e^^'' 

f[er+']-f[er] , </'.+i/2(f[e"])-4-i/2(f[e"])\ , 9K,rh:+' 



At Ax J /if At 

where the second order numerical fluxes are defined by 

A+i/2{9k) = ■"'9KMI + '"^9K,^ + ^kjl minnod (g^_, - 9K,^-l: 9K,^+l " 9K^^,9K,^+2 - 9K,^) 

with v~ = Vj if Vj < and = in other cases, while w"*" = Vj if Vj > and ti+ = if Vj is 
negative. The same numerical flux is used for 4'i+i/2{E'[g"])- Note that for simplicity, in ([5^ and 
all what follows, the discrete-velocity index j is omitted, as well as the space and time dependency 
of ly. 

Note that in ([5^ . is computed by using a discrete version of (PO)) which is explained 

below. Moreover, note that the last term of the right-hand side of (p2|) models the evolution of 
the transition function h: the new value /i""*"^ depends on the equilibrium/non equilibrium state 
of the gas in a way that will be described in section[51 In addition, we point out that when /i" = 0, 
equation ([5^ is not solved, thus the term 9xi/^2 does not lead to any computational difficulties. 
Finally, note that the stiff relaxation term of ([5^ is implicit. This allows us to use a time step 
which is independent of e. 

Now, we describe the numerical scheme for the macroscopic equation (I30|) . This equation is 
discretized according to 



g^^-gr V'.+i/2(g",gI^)-V'»-i/2(g",gI^) 

At Ax 



(33) 



where the numerical flux is an approximation of the total flux -F(g, 9k) — F{q) + {vmgK) obtained 
by the second order MUSCL extension of a Lax-Friedrichs like scheme: 

i^^-,l,2{Q\9l) = \[F{g-.gl^,) + F(gr+i, g^.+i)) - \a[g-+, gD + - af+T) (34) 
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In this relation, we set 

= {F{Q7+i.9'h+i) ± - HQ'l.gl,) T ag-) v{xT^) (35) 

where </3 is a slope limiter, a is the largest eigenvalue of the Euler system and 

n,± ^ FjQ?. 9k,) ± - F{Q?-l,9K,^-l) T «gr-i , . 

i^(er+i,g^,,+i)±aer+i-i^(er,5^^,.)Tae? ^ ^ 

where the above vectors ratios are defined componentwise. 
4.2.2 Time splitting scheme 

Here, we propose an alternative scheme based on a time splitting between the dth term and the 
other terms in the kinetic equation for gx (|3ip . We will show in the next section that this method 
preserves the positivity of the distribution function f ^ E[g] +gK under a suitable CFL condition. 

First, we solve the macroscopic equation using as in the previous scheme, where the 
numerical fluxes are defined in (p4)) . Now, for the kinetic equation on gx, the time variation of h 
only is taken into account in ([31]) to get the second step: 

1 o" 

Note that this relation can be readily simplified in 

9k7-91,^, (37) 

where, again, we point out that this equation is solved only if hf ^ 0. 

In a third step, (|3ip is discretized without the dth term, by using the same approximation as 
for the non-splitting scheme. We get: 

„n+l n+1/2 / , n+l/2x , / n+l/2x\ 

9k, ~9K,t , , „+l / 0»+l/2(gA- )~(Pi-l/2(gK ) \ _ n+1 

At [ Ax )^ 

, 0.+i/2(g[g"])-0.-i/2(g[g"]) \ 

(v — — ^ ' 

Note that for the moment, we did not mention how the new value of the transition function /i""^^ 
is defined. This is done by using some criteria that are introduced in section O Independently 
of this problem, we analyze in the following section the positivity property for the distribution 
function. 



4.3 Positivity of the distribution function for the discretized equations 

In this section, we prove that the splitting scheme p7ll38p preserves the positivity of / under a 
suitable CFL condition. Another interesting property of the model here proposed, the preservation 
of uniform flows, will be analyzed in the appendix in comparison with different coupling strategies 
proposed in the recent past [HI M EI] • 

Proposition 1. If f > and g% = h°{f - £[q°]), where q° = {mE[g]) and < /i° < 1, then 
scheme (^MSW satisfies 

fr^£[g?] + g1<,>0 
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for every n and i, provided that At satisfies the following CFL condition: 

At < — - mm —7- . (39) 

max(v,) ^,v, \h1+\gl+'/^ + £[q^]) J 

Proof. The idea is in fact to prove a stronger property: indeed, we can prove, by induction, that 
the positivity of + ^ is preserved at any time. 

First, note that this relation holds at n = 0: from the definition of 5^, we have h^£[Q^] +g% j = 

Then, we assume that this relation is satisfied for some n, and we prove that it is true for n+l. 
This is done in the following three steps. 
Step 1. 

We first use psp (where the numerical fluxes 4'i+i/2 are computed by the first order upwind scheme) 
to explicitely compute g^l and then to obtain: 

n+l I ;.n+lcr„"+ll ( I "+1/2 , l,n+lcr„nn I'^l^^ Ln+l/ n+1/2 . c\„m,\ 



V 



-At, 



Aa; 



(40) 



l + e/{iyAt) 



diti^"^ + These two terms are studied in step 2 



Now, it is clear that the sign of the left-hand side depends on the sign of g^l^"^ + h^^'^E\Q^\ and 
step 2. 

Here, we use the definition of g^+V^ (see dST])) to obtain g^t^^^^ + Zi'/^^f [£•?] = ^7^{K£[q'z]+91,i) 
which is non-negative (due to the induction assumption). Consequently, g'Z^^'^ + h^^^E[Q^] is non- 
negative. Since £[0^] > and < /i"^^ < 1, then we also have that g^^l^^ +£[Qi'] is non-negative. 

Step 3. 

Note that step 2 shows that the last three terms of the right-hand side of (|40p are non-negative. 
Consequently, (|iD|) shows that 5^+/ -|- is non- negative if At satisfies the CFL condi- 

tion (|39p . By induction, -t- g^^ ^ is non-negative for every n, and for every i and v. 

Finally, using again that £[Qi] > and < ft" < 1, we easily deduce that £[q2] +9^1 is also 
non-negative. □ 

Remark 2. If g^ i > 0, condition i39\) is less restrictive than the CFL condition At < ^^^^ ^ 
obtained with a classical semi-implicit scheme for the original BGK equation. On the contrary, 
condition h39(l becomes more restrictive if the perturbation term g^ ^ is negative. However, if we 

assume that g is small enough (e.g. £[Qf] S> 9^1^'^)) then the factor ^^^^ in \39(l is close to 1. 
Indeed: 



and i39\) reduces to the classical CFL for transport At < ^^^^^^ 
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By contrast, we justify below why we think that the non-sphtting scheme p2H33p cannot 
preserve the positivity of /. Indeed, by using similar computations as the ones we did for the 
splitting scheme, we find: 

ejv /ifAto^^„ , c-r^n n , cr^n+ii 



Now, it is clear that the sign of /"^^ depends on the signs of ft-"^^ — ft-f and g^, and hence cannot 
be controlled by a CFL condition on Ai. 

5 Localization of Fluid-Kinetic Transitions and the Dynamic 
Coupling Technique 

One of the key points in a domain decomposition technique for gas dynamics problems is to 
efficiently localize the regions where the state of the gas departs from equilibrium, so as to describe 
the solution with the appropriate microscopic model. In other words we look for an accurate 
criterion the evaluation of which is computationally inexpensive. 

Here, we propose three different criteria based on the information which can be retrieved from 
either the kinetic distribution function or from the macroscopic variables. The way the localization 
of the equilibrium and non-equilibrium regions evolves is described at the end of this section. 

5.1 Analysis of Microscopic and Macroscopic Criteria 
5.1.1 Microscopic Criteria 

In regions where the kinetic or coupled kinetic/fluid models are solved, we can use the distribution 
function to measure the fraction of gas particles which are not distributed according to a Maxwellian 
(as in [11 ). In the same way, the fractions of momentum, energy, and heat flux due to the non- 
equilibrium flux can be measured. Consequently, in every cell where /i 7^ 0, it is possible to evaluate 
the parameters 

Ai,K = {\mm. ^2,K = {v\gKiv)\), A3,K = (^|.gK(«)|), X4,k = \{v^-^9k{v))\. (42) 

Note that the first three values above will be zero if we use gx instead of |5_fs-|, while the last 
one is in general different from zero. For compatibility with the macroscopic criterion introduced 
in the sequel, the definition of A4 in (|42p has been preferred to the alternate definition X^^k = 



(\v^^gK{v)\) ■ The four parameters are computed in our code with the following quadrature 



2 

formula 



Xl,K = ^ \gK,j\^V, \2^K ^^Vj\gK,j\^V, 

j 
^a.i^^E^lff^^^l^^' X^.K^\Y.v^^gK,A^l (43) 



where gujit, x) ~ gK^x, Vj, t). In kinetic and buffer zones fli U fl^ (where h ^ 0), the discrepancy 
between the fluid and kinetic models can be measured by the following parameters 

\ n \n \n 

Pl,i,K - ' P2,i,K - ^n,,n ' P3,i,K - ^„^„ i 
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or alternatively, we can use the value of the heat flux relatively to the value of the equilibrium 
energy flux 

if 

P4,i,K = 1/ 1 I 12 E^r nl\ I ■ 

In order to define a unique variable which permits to switch from one model to the other one 
in every regime and in every region (f2i, i = 1,2,3), we choose Pi as the breakdown parameter. 
Indeed, as shown below, it is possible to estimate this quantity also in fluid regimes. By using this 
criterion, the transition function can be defined by an appropriate function that maps Pi to the 
interval [0, 1]: /i" = fiPii k)- Such a mapping is defined in section [5?2l 

5.1.2 Macroscopic criteria 

The previous analysis is quite efficient and does not induce expensive additional computations, 
since the perturbation term qk is already known in regions where the parameter P^ has to be 
computed. However, if we decide to use this criterion in the whole domain, the cost will be 
equivalent to the cost of computing the solution of the kinetic model in the whole domain. For 
this reason, it is necessary to look for others indicators that are based on the equilibrium values 
only. The most obvious one is the local Knudsen number e which is defined as the ratio of the 
mean free path of the particles Xpath to a reference length L: 

e = Xpath/ L, (46) 

where the mean free path is defined by 

kT 



X 



path 



with k the Boltzmann constant equal to 1.380062 x 10^^^ JK~^ , p the pressure and ttct^ the collision 
cross section of the molecules. The Knudsen number is determined through macroscopic quantities 
and can be computed in the whole domain with a minimum additional cost. Now, in order to take 
into account the elementary fact that, even in extremely rarefied situations, the flow can be in 
thermodynamic equilibrium, according to Bird [5], the reference length is defined as 

L = min f ^ ge \ /^^x 

\dg/dx' dgu/dx' dge/dx J 

According to [501 [IS] i the fluid model is accurate enough if the local Knudsen number is lower 
than the threshold value 0.05. It is argued that, in this way, the error between a macroscopic 
and a microscopic model is less than 5% |32j . This parameter has been extensively used in many 
works and is now considered in the rarefied gas dynamic community as an acceptable indicator. 
We notice that the local Knudsen number takes into account both the physics (with the measure 
of the mean free path and the identification of large gradients) and the numerics (through the 
approximation of derivatives on the mesh). 

In the present work we propose an alternative criterion, based on the analysis of the micro- 
macro decomposition. We will apply this criterion only in fiuid regions. For this reason, wc consider 
the equation for the non-localized perturbation g (sec (1161) ) in its discretized form 



g? , /4+i/2(.9")-</'.-i/2(g") 



At V / £ 

£[q^+'] - £[gf] ^ 0.+i/2(g[g"]) - 0.-i/2(g[g"]) 
At Ax 



(48) 
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(50) 



Let us consider a point Xi which hes in the macroscopic region at time , i.e. = 0. If, in 
addition, we assume that is close to zero in the neighboring cells (which should be true if the 
transition function is smooth enough), then we obtain 

Using this relation, we are able to evaluate the mismatch between the macroscopic fluid equations 
and the kinetic equation. In fact, note that, in the macroscopic equation (|15p we do not know how 
to evaluate the kinetic term dx < vmg >, except if we solve, at the same time, the kinetic and the 
macroscopic problem ()15II16|) . However, at point Xi, thanks to (|49|) . the perturbation term only 
depends on the Maxwellian distribution which in turn only depends on the macroscopic variables 
at the previous time step. Then, integrating (|49p over the velocity space we get: 

f vmg''+^dv = T^^^i f vm£W;+^]dv- f vmShmv] + 

^TT~7~X7Va~ <^»+i/2 ( / vvm£[Q'']dv \ - (j)i-i/2 [ [ wwm£[£>"]du 
{e/iy + At)Ax [ V^r^ / W 

where in one space dimension we have 

/ vE[Q]dv = gu, [ v^E[Q]dv = g{u^ + 39) (51) 

/ v^E[Q]dv = guiu^ + 56), [ v'^E[Q]dv ^ g{u^ + Su'^O + 59"^) (52) 

Now, the last step is to measure the ratio of the non-equilibrium fraction to the equilibrium 
one. Observe that in the one dimensional case the only non-zero moment of the non-equilibrium 
term g is the heat flux. Thus, like for the microscopic criterion (j45p . we measure the ratio of the 
heat flux to the energy flux: 

and define, as before, the value of the transition function /i" at this point as an appropriate function 
of /34: 

K = KPh). Q<K<i (54) 

In practice, we use the same function f to evaluate /i" in all the computational domain, but while 
/3Jj is defined by (|45| in kinetic regions, it is defined by ((53)) in fluid regions. 

The quantities (|45p - (|53p (which will be referred to as breakdown parameters) furnish a true 
measure of the model error, while the local Knudsen criterion is rather a physical-based criterion. 
In the numerical test section, we will compare these two strategies. 

Remark 3. In the above analysis we have discarded the convection term (vdxg)- This can be 
justified by the hypothesis of smoothly varying transitions, which means that this term is supposed 
to be small. Anyway, it is possible to take it into account. For example, through an upwind 
discretization, we obtain 



and 



vdxg"^ 



9r - 91 



Ax^ 
9^ - 9i 



- if v>Q 



V- 



- if v<0 



Ax 

Now, as before, = 0, while g][^_i or g^_i assume known values, which can be different from zero 
if the transition function h appears to be greater than zero in these cells (?i^Yi 7^ or h'^_i 0). 
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5.2 Kinetic/Fluid Coupling Algorithm 

We now describe the kinetic/fluid coupling algorithm. 

Define (3thr and P^f^^. < Pthr as the maximum errors that we can afford by using the fluid model 
instead of the kinetic one. Then: 

Assume g^, h" are known in the whole space domain at time n. 

1. Advance the macroscopic equation in time by using scheme p3|) and obtain g"^^. 

2. Compute the equilibrium parameter P2i in every space cell for which h = through relation 

(ESI. 

3. If > Pth,- then set h'^+^ = 1 which means that Xi at time t"+^ belongs to the kinetic 
region; if (3^ ^ < (3^i^^ then set h^^^ = 0, which means that Xi at time belongs to the 
fluid region. 

4- Pthr ^ P4i< Pthr then set = q'^-'Z^J!^" , which means that Xi at time belongs to 
the buffer region. 

5. Advance the kinetic equation in time by using scheme ([57|) - ([55)) and obtain g"^^ ■ 

6. Compute the equilibrium parameter /3J,- ^ in every space cell for which h =/= through relation 
(SSI). 

7. If P2,^,K > Pthr then set h'l+^ = 1 which means that Xi at time belongs to the kinetic 
region, if /?" ^ ^ < /3j^^ then set h"^^ = 0, which means that Xi at time belongs to the 
fluid region. 

8. If Pthr <Pl^K^^ Pthr thcu sct = ' ' ' 5"'' , which means that x^ at time belongs 
to the buffer region. 

9. Re-project the non equilibrium part of the distribution function through the relation (P7|) . 
Remark 4. 

• In the above algorithm the steps 2-3-4 can be substituted by equivalent steps in which the 
breakdown criterion is the local Knudsen number, with convenient threshold values. 

2. Compute the local Knudsen number ef in every space cell for which h = through 
relation |^6p . 

3. Ife'2 > etrh then set K^+^ = 1 which means that Xi at time belongs to the kinetic 
region, if < e^^^ then set h"^^ — 0, which means that Xi at time t"^^ belongs to the 
fluid region. 

4- If ^trh — — ^trh then set h^l^^ — ^^'^J^^';'' , which means that Xi at time t"^^^ belongs 
to the buffer region. 

In the next section, we report comparisons between using the Knudsen number and the new 
breakdown parameter. 

• At the beginning of the computation the full domain is supposed in thermodynamical equilib- 
rium. During the computation, kinetic regions are created. Some of these regions can become 
even one cell thick, merge or split. The transition function can also pass from to 1 and 
vice-versa in a single time step and with jumps in space. Every step is completely automatic 
in each .simulation, no additional parameters are used. 
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• Compared to the previous work JTTI in which the buffer regions where fixed once for all at 
the beginning of the computation, the present strategy consists in making the buffer regions 
dependent on the current state of the gas through the functions '^^'^ thresholds 

values. This modification leads to a considerable improvement in accuracy, flexibility and 
usability of the method. 

6 Numerical tests 
6.1 General setting 

In this section, we present several numerical results to highlight the performances of the method. 
By using unsteady test problems, we emphasize the deficiencies of the static decomposition method. 
As in [11] we start with an unsteady shock test problem. Even in this simple situation, a standard 
static domain decomposition fails in its scope. Indeed, the shock moves in time. Thus in rarefied 
regimes, it is necessary to use a kinetic solver in the full domain, which turns to be a quite inefhcient 
method. On the other hand, with our algorithm, we reduce the computationally expensive regions 
to a small zone compared to the full domain. 

Next, we use our scheme to compute the solution of the Sod test. Here some new difficulties 
arise. Indeed, contact discontinuities and rarefaction waves appear but, as described below, the 
method efficiently deals with such situations. 

Finally, in the third test, a blast wave simulation is performed. In spite of the complexity of the 
solution, the algorithm shows a very good behavior, creating, deleting or merging zones together 
and obtaining fast and precise results. 

In order to obtain the correct equation of state with only one velocity-space dimension, we use 
the following model: 

^* ( G ) + ( G ) = ^ ( tLp^G ) ■ 

It is obtained from the full three-dimensional Boltzmann-BGK system by means of a reduction 
technique [16j . In this model, the fluid energy is given by 

E = [ {^v^F + G)dv. 

This model permits to recover the correct hydrodynamic limit given by the standard Euler system 
even with a one-dimensional velocity space. 

The collision frequency is given by i' = t^^ = (^)^^ where fi — C-O'^. We choose gas hydrogen 
for our simulations. Thus C = 1.99 x 10"^^ lo = 0.81 and R = 208.24 [2J. 

In all tests the time step is given by the minimum of the maximum time steps allowed by the 
kinetic and fluid schemes. This means that no attempts have been made to try to increase the 
computational efficiency by means of a reduction of the number of necessary effective steps for 
the less restrictive scheme, by, e.g. freezing one model in time, while the other one is advanced. 
Indeed, such a reduction still requires further investigations. Thus, the global speed-up is only due 
to the reduction of the sizes of the kinetic and buffer regions inside the domain. This reduction 
is achieved through a correct prediction of the evolution of the transition function and the use of 
efficient criteria for the determination of the equilibrium regions. We point out that no a-priori 
choices on the dimension of the buffer and position of the different regions are done in all the tests. 
Instead all the procedure is automatic and determined by the step by step algorithm presented in 
the previous section. For all tests, we set f3thr = 10~^ and (i^f^^ = 10~^. 
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6.2 Unsteady Shock Tests 



We consider an unsteady shock that propagates from left to right in the computational domain 
X = —20 m, X — 2Q m. The shock is produced pushing the hydrogen gas against a wall which 
is located on the left boundary. We consider that the particles are specularly reflected and that 
the wall instantaneously adopts the temperature of the gas. This effect is numerically simulated 
by setting macroscopic variables in ghost cells (two cells for a second order scheme) beyond the 
left boundary with parameters g, T equal to the values of the first cell while the momentum is 
set opposite. In the non equilibrium case gx is also different from zero in the ghost cell, and is 
equal to a specularly reflected copy of in the first cell. At the right boundary, we mimic the 
introduction of the gas by adding two ghost cells where, at each time step, the initial values for 
density, momentum and energy are fixed. The gas is supposed in thermodynamic equilibrium, 
which implies that gx = 0. The computation is stopped at the final time t = 0.04 s. There are 
300 cells in physical space and 40 cells in velocity space. We do not use a finer mesh because the 
scheme is second order. Symmetric artificial boundaries in velocity space are fixed at the beginning 
of the computation through the relation Vb — ±Ci max('\/RTvy), where R is the gas constant, C'l 
is a parameter normally fixed equal to 4 and Tw is a temperature set equal to the maximum 
attainable temperature, which is obtained by supposing that all the kinetic energy is transformed 
into thermal energy. The transition function h is initialized as h — (fluid region) everywhere. 

In the first test the initial conditions are g — 10^^ kg/m^ for the mass density, u = —900 
m/s for the mean velocity and T = 273 K for the temperature. In Figure [1] we have reported the 
mass density on the left and the velocity on the right. ^^From top to bottom, time increases from 
t = 10"^ s (top) to i = 4 X 10~2 s (bottom), with t ^2x 10"^ s middle top and i = 3 x 10"^ s 
middle bottom. In Figure[5]we have reported the temperature on the left, the transition function, 
the heat flux and the local Knudsen number on the right. From top to bottom the same instants 
of time as for the previous Figure are shown. In the plots regarding the macroscopic variables we 
reported the solution computed with our algorithm (mic-mac in the legend of the Figures), the 
solution computed with a kinetic solver and the solution computed with a macroscopic fluid solver. 
Magnifications of the solutions close to non equilibrium regions are given for clarity. 

As soon as the simulation starts on the left boundary, the transition function h increases from 
zero to one, which means that the solution is computed with the kinetic scheme, while in the rest 
of the domain the solution is still computed with the fluid scheme {h = 0). When the shock starts 
to move towards the right, we notice a splitting of the kinetic region. One very narrow region still 
continues to follow the shock and one remains close to the left boundary. Once that the threshold 
values of the breakdown parameters /? and Pk are fixed, the procedure automatically determines 
the sizes of the kinetic and buffer regions. 

We repeat the simulation with a lower initial density g = 10~^ kg/m^. This yields different 
results which are reported in Figure [3] for the density and velocity and in Figure [4] for the temper- 
ature, local Knudsen, heat flux and transition function. The results are obtained with the same 
criteria as in the previous test case. Again at the beginning h is set equal to zero (fluid), but now 
the shock is much less sharp and the non equilibrium region becomes larger. 

We observe that in this first test the discrepancy between the fully kinetic solver and the 
coupling strategy is not perceivable, while the computational time is reduced in proportion to 
the ratio between the areas of the kinetic and buffer regions compared to the entire domain. 
Thus, in the first case we have a reduction of approximately 70% of the computational time while 
in the second one the reduction is only the 20%. Further reductions are possible through code 
optimization. 
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6.3 Sod Tests 



In these second series of tests, we consider the classical Sod initial data in a domain which ranges 
from —20 m to 20 m. The numerical parameters are the same as for the previous examples. 
Thus respectively 300 and 40 mesh points in physical and velocity space are used. Symmetric 
artificial boundary condition are fixed in velocity space at the same points Uf, — ±Ci max(^/'RTvv'), 
while Neumann boundary condition are chosen both for the kinetic (if necessary) and fluid models. 
Finally the simulations are initialized with a thermodynamic equilibrium with h = Q and qk = 
everywhere. 

In the first case, we take the following initial conditions: mass density ql — 2 y. 10^^ kg/m'^, 
mean velocity ul — ^ m/s and temperature Tl — 273.15 K if —20 < x < m, while qr — 
0.25 X lO"'^ kg/m^, ur = {) m/s, Tr = 218.4 K if < x < 20 m. The results are reported in Figure 
[5] for the density (left) and velocity (right) and in Figure [5] for the temperature (left), Knudsen 
number, heat flux and transition function (right). Snapshot at increasing times are displayed top 
to down, corresponding successively to i = 6 x lO^'^ s, then i = 1.2 x 10^^ s, t = 1.8 x 10^^ s and 
finally t — 2.4 x 10~^ s. Again we provide magnifications of the solution close to non equilibrium 
zones in order to highlight the differences between the three different schemes: the macroscopic 
and kinetic ones and the coupling strategy (mic-mac in the legend). We observe that due to the 
initial shock, a kinetic region appears immediately and starts to grow in time, but as soon as the 
different non equilibrium regions separate, the kinetic region itself splits into three: one around 
the rarefaction wave, one around the contact discontinuity, and one around the shock. Even with 
magnifications it is not possible to perceive differences between the kinetic model and the coupling 
strategy, even though the kinetic regions remain very tiny, permitting a fast computation. 

The simulation is repeated, with a lower initial density = 5 x 10^^ kg/m^ and qr — 
0.75 X 10~^ kg/m'^, and the results are displayed in Figure [7] for the density and velocity and 
in Figure [8] for the temperature, heat flux, local Knudsen and transition function. The same 
qualitative features as in the previous test can be observed, the only difference being that the 
kinetic regions are thicker, which means that the non equilibrium zone is larger. This is clearly 
visible from the plots of the macroscopic quantities: the difference between the kinetic and fluid 
models is significant in a large portion of the domain. 

We finally observe that compared to |11| , in which a similar scheme was developed, we are able to 
capture small discrepancies between the kinetic and macroscopic models in very tiny regions. This 
turns to be a much more efficient use of the domain decomposition technique. The computational 
time reduction is of the order of 70% and 60% respectively for the two tests compared to a kinetic 
solver. 

6.4 Blast Wave Tests 

In this paragraph we present two interacting symmetric blast waves in hydrogen gas. The domain 
ranges from x = Omtoa; = lm, while the numerical parameters in terms of mesh and velocity 
space boundaries are the same as in the previous tests. The physical boundaries are represented 
by two specularly reflecting walls, on which impinging particles are re-emitted in the opposite 
direction with the same velocity (in magnitude). Mass is conserved at the walls which additionally 
are supposed to adopt the gas temperature instantaneously. These effects are obtained like in 
the unsteady shock test with two ghost cells (four for a second order scheme), in which the same 
macroscopic values as those of the first and last cells are imposed, except for momentum which 
changes sign. The perturbation function can in general be different from zero and assumes the 
same values of its corresponding counterpart in the first and last cell, with a sign change in the 
velocity variable. At the beginning, we set h = and = everywhere. 
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In the first test, the initial data are: 

Q = IQ-^kg/m^ u = 200 m/s T = 10000 K if x < 0.1 

g = IQ-^kg/m^ u = -200 m/s T = 10000 K if x > 0.9 

g^lO~^kg/m^ u = Om/s T = 50 K if OA < x < 0.9 

The results in terms of the density and mean velocity are reported in Figure [9l while the temper- 
ature, local Knudsen, heat flux and transition function are reported in Figure [101 The displayed 
results are for increasing times t = 10"'' s to t — 4 x 10^"* s from top to bottom. Again we 
plot the kinetic, the fluid schemes and the coupling strategy (micmac scheme) in each Figure and 
magnifications close to non equilibrium regions are provided. 
In the second test the initial density is decreased, and we use: 

Q = lO^^kg/m^ u = 200 m/s T = 10000 K if x < 0.1 

g = lO^^kg/m^ u = -200 m/s T = 10000 K if x > 0.9 

g=10^'^kg/m^ u^Om/s T ^ 50 K if 0.1 < x < 0.9 

The results obtained with this second set of data are reported in Figure [Tl] and [121 We observe 
that starting from a situation where the fluid model is used almost everywhere we end up in the 
opposite situation where the kinetic model is used in the whole domain(/i = 1 Va;). Thus, while 
a static domain decomposition technique leads to similar computational times as a fully kinetic 
resolution, the coupling strategy leads to a speed up of about 40% for the first test and 30% for 
the second test, compared to a fully kinetic solver. 

7 Conclusion 

In this paper we have presented a moving domain decomposition method which provides an efficient 
way to deal with multiscale fluid dynamic problems. Regions far from thermodynamical equilibrium 
are treated with a kinetic solver. The method is based on the micro-macro decomposition technique 
developed by Degond-Liu-Mieussens [lOj in which macroscopic fluid equations are coupled with 
a kinetic equation which describes the time evolution of the perturbation from equilibrium. The 
method consists in splitting the distribution function into an equilibrium part and a non-equilibrium 
part, together with the introduction of buffer zones and transition functions as proposed in [10], 
[8] and [11] to smoothly pass from the macroscopic model to the kinetic model and vice versa. 

In order to build up an efficient method that can be used in a wide spectrum of situations, 
and by contrast to [10], we consider the possibility of moving the different domains like in [TT] . 
using however, a decomposition technique which shows enhanced performances. Moreover, we have 
developed a scheme which is able to automatically create, cancel and move as many kinetic, fluid 
or buffer regions as necessary. The method relies on the proper combination of the two equilibrium 
criteria we have identifled and on a priori tolerance that we decide to accept. An important point 
also resides in the introduction of a new criterion for the breakdown of the fluid model, which is 
able to measure the discrepancy between the kinetic and the fluid model in a much more accurate 
way then the mere Knudsen number. Finally, we have proved that the coupling strategy preserves 
positivity under a CFL condition, and the uniform flows, also in the fully discrete case. 

The last part of the work is devoted to several numerical tests. The results clearly demonstrate 
the advantages of this method over existing ones. The method captures the correct kinetic behav- 
iors even in transition regions and provides signiflcant improvements in terms of computational 
speedup while maintaining the accuracy of a kinetic solver. 
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In the future we will extend the method to make it consistent with the Navier-Stokes equations 
instead of the Euler model. This step can further improve the technique allowing very narrow 
kinetic zones and providing considerable speed-up while maintaining accuracy. We will also explore 
the use of Monte Carlo techniques for the full Boltzmann equation and that of time sub-cycling 
for the two models. We finally observe that the computational speed-up will significantly increase 
for two or three dimensional simulations, which we intend to carry out in the future. To conclude, 
because multiscale effects are very important also in many others fields we plan to extend our 
method to other models. 



A Preservation of uniform flows 

Preservation of uniform flows is a very important property, which prevents oscillations to appear 
when the transition region is located in a domain where the flow is smooth. In this appendix we 
compare the properties of the present micro-macro coupling strategy to those of the decomposition 
methods of [8l[9l|TT] regarding preservation of uniform flows. 

In TO] , it has been demonstrated that the micro-macro model is able to preserve uniform flows 
in the continuous case. This property is also true for the decomposition used in [81 ITT], but 
only when the collision operator has specific properties (which are satisfied by the Boltzmann and 
BGK operator) . In this appendix, we show that the present micro-macro coupling strategy is able 
to preserve uniform flows also in the discrete case independently of the choice of the discretization 
scheme. We observe that this property does not hold in the general case for the decompositions 
used in [HI [9l [11] . The satisfaction of this property by the present micro- macro coupling strategy 
constitutes a very big advantadge of this method over the previous ones [H [9l [11] . 

As an example, we consider the decomposition used in fTD , which reads 

^ + (1 - h)d,F{Q^) + {1- h)d^{vmfR) - -Qdth (55) 

dtjR + hvd^R + hvd,E[Q] = h-^{E[Q\ -f) + fdth, (56) 

where the distribution function is defined hy f = Jr + E[Q]^], Jr — hf and E[qi^] = (1 — h)E[Q]. 
In this model the solution of the full kinetic problem is given by /rUxG^i, by E[qi] ii x E il2 
and by E[qi^] + fRif x G Os. This is due to the fact that fu = for x G fl2, E[qi^] = for x G ili, 
while in fis they are both different from zero and the global solution is obtained as a sum of the 
two partial solutions. We refer to the above cited papers for details. 

Here we recall that, in |SI[H1[II], small oscillations appear inside the transition regions except 
when the scheme used for the fluid part is an exact discrete velocity integration of the scheme used 
for the kinetic part (in this case, we say that the two schemes are 'compatible', and are 'incompat- 
ible' otherwise). These oscillations appear even in situations where preservation of uniform flows 
is true for the continuous model. To circumvent this problem, in [H [9l [H], we used a standard 
shock-capturing scheme (such as e.g. the Godunov scheme) for the Euler equations in the pure 
fluid region (i.e. h = 0), but we converted to a compatible scheme with the discretization of the 
kinetic model inside the buffer zones (see [H] for details). However, this choice introduces some 
implementation difficulties and reduces the performances. Indeed, a compatible scheme with the 
discretization of the kinetic model has intrinsically the same cost as the full kinetic solver, and the 
coupling strategy is twice more costly than the mere kinetic model in all the buffer region. 

In order to prove the property that uniform flows are not preserved by the decomposition (j55I 
[56]) in the discrete case, we focus on the first one of the two equations. The same considerations 
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hold for the other one. If the initial data is such that / = E[q] we have 
+ (1 - h)d^F{QL) + (1 - h)d^{vmfR) + Qdth = 



dQL 



dt 

= (1 - /i) (^ll + (1 - h)d^F{Q) + hd^{vmE[Q\) - h'F{Q) + h'{vmE[g]) 
= (1 - M + (1 - h)d^F{Q) + hd4vmE[Q]) 

In the above equation, the time derivative with respect to h disappears and so does the flux, using 
that F{q) = {vmE[Q]). In the continuous case it is also true that 

(1 - h)d^F{Q) + hd^{vmE[Q]) = d^F{Q) = d^{vmE[g]) (57) 

and so, uniform flows are preserved. However, when we discretize the fluxes according to dxF{Q) — 
('/'z+i/2(e)-0»-i/2(£'))/Ax a,nddx {vm£{Q)) = (-0j+i/2(^[e])-V',-i/2(^[e]))/Aa:, the equality ^ 
does not hold anymore. In fact, in the general case we have 

fl_h\( '^»+i/2(g) - 0»-l/2(g) \ ^ ^ / V'»+i/2(g[g]) - V'»-l/2(g[g]) \ ^ 
\ Ax J \ Ax J 



Ax J \ Aa; 

Thus, if two incompatible numerical schemes are used to discretize the kinetic and fluid fluxes, oscil- 
lations in the solution can appear as documented in [llj . However, we observe that, using the same 
numerical flux is not sufficient to ensure preservation of uniform flows through the decomposition 
(|55ll56p . To that aim, consider a generic discretization of the coupled system (|55ll56p : 

- (1 " ^^)|^ E (V'^+l/2(/fc,fl) - ^^~l/2{flR)) A^, (58) 

k 

fk^R ^ fk,i,R - (4+1/2 (/fc,_R.) - (t'i-l/2{fkM)) 

-h~ (4+1/2 (ffc[e£]) - 4-1/2 (ffe[e2])) (59) 

e 

where a discrete velocity model has been used to resolve the kinetic equation ([59]) with m.k the 
discretized collision invariants. The function 4±i/2j V'i±i/2 are two different generic numerical 
fluxes while, for simplicity, the function h is considered constant in time. The initial data are 
q1 = /o = E[eol fl, = hjo, qL, - (1 - h,)Q^ and E\qI^^ - (1 - h;)E\Q^\ Vz. Now, 
supposing the flow uniform at time n we will see that not every numerical scheme ensures a 
uniform flow at time n + 1. To this aim, if we integrate equation (|59p multiplied by the collision 
invariants over the velocity space, we can rewrite the coupled numerical schemes (|58ll59p as 
follows: 

qTl' ^Y.^3QIl. (60) 
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i±Il K 



(61) 



with — X]/c mfc/ffl^^^i I ^-nd /I the length of the stencils in physical space and K in velocity 
space. The symbols Aj and Bj^k represent the weights determined by the particular choice of the 
numerical schemes. Without loss of generality, suppose that I = II. Then, we have that 



n+l 
tii.L 



E 



E 



) k 



E^.^" + E 



E 



Y^B,.^mu!lA-A,Ql^^ 




i-R 



(62) 



which means that we do not have preservation of uniform flows except in some particular cases, 
such as, for instance, when the numerical schemes used to discretize the two equations (|58ll59p are 
compatible. Therefore, such compatible schemes are needed in all buffer zones to make sure that 
oscillations in the solutions will be avoided. 

Instead if we consider the present micro-macro coupling strategy with the same initial data 
/ — E\q\-, the following property holds: 

Proposition 2. // the initial condition f^>Oisa constant equilibrium E[g^], then Q = and 
9k — h{f — E[q\) = are solutions of the micro-macro model h2S^24^ , and E[q\ + gx — E[q'^]- 
In other words, the kinetic/fluid solution of the micro-macro model is exactly the solution of the 
original kinetic model. 

Indeed, for the micro-macro decomposition, the total flux is independent of h in equilibrium 
regimes and is not obtained as a sum of two complementary terms weighted by the function h. 
It follows directly that the micro-macro method preserves uniform flows even in the discrete case 
independently of the choice of the numerical scheme. 
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Figure 1: Unsteady Shock 1: Solution at t = 1 x 10^^ top, t = 2 x 10^^ middle top, t 
3 X 10^^ middle bottom, t — 4 x 10^^ bottom, density left, velocity right. The small panels are 
magnification of the solution close to the shock. 
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Tamperalure Unsteady Shock TbsI [or t-0.01 s 
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Figure 2: Unsteady Shock 1: Solution at t = 1 x 10"^ top, i = 2 x 10"^ middle top, i = 3 x 10"^ 
middle bottom, t — A x 10~^ bottom, temperature left, transition function, Knudsen number and 
heat flux right. The small panels are a magnification of the solution close to the shock. 
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Figure 4: Unsteady Shock 2: Solution at < = 1 x 10"^ top, t = 2x 10"^ middle top, t = 3x 10"^ 
middle bottom, t = A x 10^^ bottom, temperature left, transition function, Knudsen number and 
heat flux right. 
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VslDclty Sod Test tor UO.D06 s 




Figures: Sod Test 1: Solution at < = 0.6 x lO^Hop, < = 1.2 x 10"^ middle top, i = 1.8x 10"^ j^iddle 
bottom, t — 2.4 X 10^^ bottom, density left, velocity right. The small panels are a magnification 
of the solution close to non equilibrium regions. 
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Temperature Sod Test tor (^0.006 s 
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Figure 6: Sod Test 1: Solution at t ^ 0.6 x 10"^ ^^p^ t = 1.2 x 10"^ middle top, i = 1.8 x 10"^ 
middle bottom, t = 2.4 x 10^^ bottom, temperature left, transition function , Knudsen number 
and heat flux right. The small panels are a magnification of the solution close to non equilibrium 
regions. 
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□snsily Sod Test tor 1^ 




-20 -15 -10 



10 15 20 -20 -15 -10 -5 



15 20 



Figure?: Sod Test 2: Solution at < = 0.6 x lO^Hop, < = 1.2 x 10"^ middle top, i = 1.8x 10"^ j^iddle 
bottom, t — 2.4 X 10^^ bottom, density left, velocity right. The small panels are a magnification 
of the solution close to non equilibrium regions. 
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Temperature Sod Test tor (^0.006 




Figure 8: Sod Test 2: Solution at i = 0.6 x 10"^ ^p, t = 1.2 x 10"^ middle top, i = 1.8 x 10"^ 
middle bottom, t = 2.4 x 10^^ bottom, temperature left, transition function, Knudsen number 
and heat flux right. The small panels are a magnification of the solution close to non equilibrium 
regions. 
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Density Blast Wave Test lor t=1e^E 



Velocity Blast Wave Test lor t=1e'' s 
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Figure 10: Blast Wave Test 1: Solution at i = 1 x 10"'' top, t = 2x 10'^ middle top, t = 3 x 10^"* 
middle bottom, i = 4 x lO"'* bottom, temperature left, transition function, Knudsen number and 
heat fl-ux right. The small panels are a magnification of the solution close to non equilibrium 
regions. 
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Density Blasl Wave Test for 1=1 e"'s 
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Figure 12: Blast Wave Test 2: Solution at i = 1 x 10"'' top, t = 2x 10'^ middle top, t = 3 x lO^"* 
middle bottom, t — 4 x 10"'* bottom, temperature left, transition function, Knudsen number and 
heat flux. 
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